// --------------------------------------------------------------------
//
//  File:        TAk12D3.rtf [Twisted Ate pairing for twisted Edwards curves with sextic twists]
//  Date:        Oct. 10, 2018
//  Last update: Oct. 10, 2018
//  Description: Magma implementation of the twisted ate pairing on twisted Edwards curves with sextic twists.
//               Embedding degree k = 12, CM discriminant D = 3
//
//  (c) 2018, Georgios Fotiadis & Chloe Martindale
//                 
//  Georgios Fotiadis, University of the Aegean, Greece, email: gfotiadis@aegean.gr
//  Chloe Martindale, Technische Universiteit Eindhoven, the Netherlands, email: chloemartindale@gmail.com
// --------------------------------------------------------------------
//  Polynomial Family
//  p(x) = (x^6 − 2*x^5 + 2*x^3 + x + 1)/3
//  t(x) = x + 1
//  r(x) = x^4 - x^2 + 1
// --------------------------------------------------------------------
//  Twisted ate pairing: a(P,Q) = [ f_{T2,P}(Q) ]^((p^k - 1)/q)
// --------------------------------------------------------------------

//Weierstrass to Montgomery
MontPoint:=function(P, alpha, s)
  x:=s*(P[1] - alpha);
  y:=s*P[2];
  MP:=Vector(2,[x,y]);
  return MP;
end function;

//Montgomery to Edwards
EdPoint:=function(P)
  X:=(P[1]/P[2]);
  Y:=(P[1] - 1)/(P[1] + 1);
  Z:= 1 ;
  W:=(X*Y)/Z;
  EP:=Vector(4,[X,Y,W,Z]);
  return EP;
end function;

//Miller's function point doubling: h_{P,P}(Q) 
hDBL:=function(P, Eda, Edd, eta, theta)
  X1:=P[1]; Y1:=P[2]; W1:=P[3]; Z1:=P[4];
  A:=X1^2;
  B:=Y1^2;
  C:=Z1^2;
  D:=Eda*A;
  E:=B + D;
  F:=2*C - E;
  G:=(X1 + Y1)^2 - A - B;
  H:=(Y1 + Z1)^2 - B - C;
  X3:=G*F;
  Y3:=E*(B - D);
  Z3:=E*F;
  W3:=G*(B - D);
  CX:=H - 2*D;
  CY:=(X1 + Z1)^2 - A - C - G;
  CW:=Edd*((X1 + W1)^2 - A) - C - E;
  R:=Vector(4,[X3,Y3,W3,Z3]); 
  h:=CX*eta + CW*theta + CY;
  return R, h;
end function;

//Miller's function point addition h_{P,S}(Q) 
hADD:=function(P, S, Eda, Edd, eta, theta)
  X1:=P[1]; Y1:=P[2]; W1:=P[3]; Z1:=P[4];
  X2:=S[1]; Y2:=S[2]; W2:=S[3]; Z2:=S[4];
  A:=X1*X2;
  B:=Y1*Y2;
  C:=Z1*W2;
  D:=Z2*W1;
  E:=W1*W2;
  F:=(X1 - Y1)*(X2 + Y2) - A + B;
  G:=B + Eda*A;
  H:=D - C;
  I:=D + C;
  X3:=I*F;
  Y3:=G*H;
  Z3:=F*G;
  W3:=I*H;
  CX:=(W1 - Y1)*(W2 + Y2) - E + B + H;
  CW:=X2*Z1 - X1*Z2 - F;
  CY:=(X1 - W1)*(X2 + W2) - A + E;
  R:=Vector(4,[X3,Y3,W3,Z3]); 
  h:=CX*eta + CW*theta + CY;
  return R, h;
end function;

//final exponentiation k = 12, D = 3
finalexp:=function(f, x0, p, q, Fp)
  t0:=1;
  f:=Frobenius(f, Fp, 6)/f; 
  f:=Frobenius(f, Fp, 2)*f; 
  x1:=x0^2;
  l:=((x1 - 2*x0 + 1) div 3);
  y:=f^l;
  t:=Frobenius(y, Fp, 3);
  t0:=t0*t;
  t:=Frobenius(y^x0, Fp, 2);
  t0:=t0*t;
  y:=y^(x1 - 1);
  t:=Frobenius(y, Fp, 1);
  y:=y^x0;
  t0:=t0*t*y*f;
  return t0;
end function; 

//Twisted ate pairing calculation
TwistedAatePairing:=function(P, Te, q, Eda, Edd, p, k, z0, eta, theta, Fp)
  f:=1;
  T:=P;
  i:=Floor(Log(2,Te)) - 1; 
  s:=Intseq(Te,2); 
  while i ge 0 do
    T,h:=hDBL(T, Eda, Edd, eta, theta); 
    f:=h*f^2; 
    if s[i+1] eq 1 then
      T,h:=hADD(P, T, Eda, Edd, eta, theta); 
      f:=h*f; 
    end if;
    i:=i - 1;
  end while;
  f:=finalexp(f, z0, p, q, Fp);
  return f;
end function;

//MAIN PROGRAM
t2:=0; LOOP:= 100 ;
D := 3 ; k := 12 ;                 
p := 13134002065465082451978965994689008063352936155273231976285809101659121364125466002052917860807951110122803152950817 ;              
q := 115792089237317700346157840501962347976585376275497230292808382329089323460673 ;      
h := 113427455640313558237160070382243949568 ;    
z0:= 18446744073709611553 ; 
a:= 0 ; 
b:= 1 ;

s:=k div 6;
Fp := GF(p);
P<x> := PolynomialRing(Fp);

E:=EllipticCurve([Fp!(a), Fp!(b)]); 
f:=x^3 + a*x + b; 
n:=h*q; 
t:=p + 1 - n; 
T:= t - 1; 
Te:=T^2 mod q;

ROOTS:=Roots(f);
ALPHA:=Vector(3,[ ROOTS[1][1],ROOTS[2][1],ROOTS[3][1] ]);

for i:= 1 to 3 do
  test,sqrt:=IsSquare(1/(3*ALPHA[i]^2 + a));
  if (test eq true) then
    ss:=sqrt;
    alpha:=ALPHA[i];  
    break i;
  end if;
end for;

//Montgomery Curve
A:=Fp!(3*alpha*ss); 
B:=Fp!(ss);
//Twisted Edwards Curve
Eda:=Fp!((A + 2)/B); 
Edd:=Fp!((A - 2)/B); 

M:=Fp!(2*(Eda + Edd)/(Eda - Edd)); 
N:=Fp!(4/(Eda - Edd)); 

P<z> := PolynomialRing(Fp);
g:=IrreduciblePolynomial(Fp,s); 
Fps<u> := ExtensionField<Fp,z|g>;
R<zz> := PolynomialRing(Fps);

for i1:= 1 to 5 do
  for i0:= 0 to 5 do
    fX:=i1*u + i0; 
    FX:=zz^(k div s) - fX;
    if IsIrreducible(FX) eq true then
      Fpk<v> := ExtensionField<Fps,zz|FX>; 
      x0:=Roots(FX, Fpk); 
      w:=x0[1][1];
      break i1;
    end if;
  end for;
end for;

W:=w; 
bt:=(-M^3*N^3)/(27*W^6);
Et := EllipticCurve([Fps!(0),Fps!(bt)]);
ht:=#Et div q; 
if IsZero(#Et mod q) eq false then
  W:=w^5; 
  bt:=(-M^3*N^3)/(27*W^6);
  Et := EllipticCurve([Fps!(0),Fps!(bt)]); 
  ht:=#Et div q; 
end if;

bk:=Fpk!((-M^3*N^3)/(27));
Ek:=EllipticCurve([Fpk!(a), Fpk!(bk)]); 
Zk<xx> := PolynomialRing(Integers());
NN := Resultant(xx^k - 1, xx^2 - t*xx + p);
hk := NN div q^2;

for i:= 1 to LOOP do
//create two points PW and Pa in E(Fp) of order q (two points because of the bilinearity check)
G:=Random(E); 
nP:=h; 
P:=nP*G; 
Pa:=P;                               
nP:=Random(1, q - 1); 
PW:=nP*P;   

MP:=MontPoint(PW, alpha, ss);
EP:=EdPoint(MP);
MPa:=MontPoint(Pa, alpha, ss);
EPa:=EdPoint(MPa);

//create two points QW and Qa in Et(Fp^2) of order q (two points because of the bilinearity check)
G:=Random(Et); nQ:= ht ; 
Q:=nQ*G;  
QW:=Q; 
xQ:=QW[1]*W^(2); 
yQ:=QW[2]*W^(3); 
E1:=Ek![xQ,yQ]; 

XQ:=Frobenius(xQ, Fp, 1); 
YQ:=Frobenius(yQ, Fp, 1); 
D1:=Ek![XQ,YQ]; 
QW:=D1-E1; 

xQ:=QW[1]*W^(-2);
yQ:=QW[2]*W^(-3);
QW:=Et![xQ,yQ]; 

nQ:=Random(1, q - 1); 
Qa:=QW; 
QW:=nQ*QW; 

eta:=N*(3*QW[1]*W^2 + 3*N - M*N)/(6*W^3*QW[2]);
theta:=N*(3*QW[1]*W^2 - 3*N - M*N)/(6*W^3*QW[2]);
etaa:=N*(3*Qa[1]*W^2 + 3*N - M*N)/(6*W^3*Qa[2]);
thetaa:=N*(3*Qa[1]*W^2 - 3*N - M*N)/(6*W^3*Qa[2]);
nPQ:=nP*nQ;

t1 := Cputime();
twistedate:=TwistedAatePairing(EP, Te, q, Eda, Edd, p, k, z0, eta, theta, Fp);
t2:=t2 + Cputime(t1); 
end for;
t2:=t2/LOOP;
printf"average time for %o pairings = %o\n", LOOP, t2;

//Bilinearity check
TwistedAatePairing(EP, Te, q, Eda, Edd, p, k, z0, eta, theta, Fp) eq TwistedAatePairing(EPa, Te, q, Eda, Edd, p, k, z0, etaa, thetaa, Fp)^(nPQ);